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ABSTRACT 

We use high resolution numerical simulations to investigate the adiabatic contraction of dark matter 
halos with a Hernquist density profile. We test the response of the halos to the growth of additional 
axisymmetric disk potentials with various central concentrations and the spherically symmetric potential 
of a softened point mass. Adding the potentials on timescales that are long compared to the dynamical 
time scale of the halo, the contracted halos have density profiles that are in excellent agreement with 
analytical predictions based on the conservation of the adiabatic invariant M{r)r. This is surprising 
as this quantity is strictly conserved only for particles on circular orbits and in spherically symmetric 
potentials. If the same potentials are added on timescales that are short compared to the dynamical 
timescale, the result depends strongly on the adopted potential. The adiabatic approximation still works 
for disk potentials. It does, however, fail for the central potential. 

Subject headings: dark matter — galaxies: formation — galaxies: halos 



1. INTRODUCTION 

The hierarchical clustering model is the paradigm for 
galaxy formation today. In this cosmological scenario, 
structure forms around peaks of primordial dark matter 
density fluctuations. The baryonic matter, which can dis- 
sipate energy through radiation, cools and falls into the 
center of its surrounding dark halo. The question of how a 
spherical mass distribution, e.g. a galactic bulge, will re- 
spond to the growth of mass in its center has been adressed 
by Barnes & White (1984), hereafter BW84. They de- 
vised a simple recipe for predicting the density profile of a 
contracted spherical density distribution. This recipe was 
used by Blumenthal et al. (1986) to examine contracted 
dark halos, assuming that the baryonic disk forms in the 
center so slowly that the time for the increase of mass in- 
side an orbit of a dark particle is long compared to its 
orbital period. In a slowly varying potential the action in- 
tegral ji = J Pidqi is a conserved property of the particle 
orbit, called an adiabatic invariant (Binney & Tremaine 
1987). Here, j is the action, p and q are the phase space 
coordinates of the dark particle. As a first approxima- 
tion, BW84 assumed a spherical density distribution with 
particles moving on circular orbits. In this case, the ra- 
dial action integral simplifies to the conservation of angu- 
lar momentum. With L = mvr and the circular velocity 
v'^{r) = GM{r)/r we get the adiabatic invariant M{r)r. 
Given that we know the final baryonic matter distribution, 
e.g. an exponential disk-like profile Mb{r) and the initial 
dark halo distribution Mi{r), we can construct the final 
dark matter distribution Mf{r) according to 

rf[M,{rf) + Mf{r)]=rNh{r), (1) 
where rf'is the final radius of a dark particle. This approx- 
imation is strictly valid only if the initial mass distribution 
is spherically symmetric. The mapping between initial and 
final radius is unique, because the dark particles are dis- 
sipationless and their circular orbits do not cross. The 
adiabatic compression has been used widely in estimating 



rotation curves in semi-analytical galaxy models (Ryden 
& Gunn 1987; Ryden 1988, 1991; Flores et al. 1993; Mo, 
Mao, & White 1998), in investigating the origin of the 
TuUy-Fisher relation (Courteau & Rix 1999), in analysing 
the core structure of dark matter halos (van den Bosch 
& Swaters 2001; Marchesini et al. 2002) and in the for- 
mation of gaseous disks in cosmological N-body halos (M. 
Steinmetz 2001, private communication). The adiabatic 
contraction assumption has been tested qualitatively in a 
low-resolution study by Blumenthal et al. (1986). 

It is surprising that equation (1) should hold in realistic 
situations, where a non-spherical galactic disk is added to 
a halo of particles which move on eccentric or chaotic or- 
bits (Valluri & Merritt 2000). In this paper we investigate 
the validity of equation (1), using high resolution numeri- 
cal simulations. Adding a disk potential to an equilibrium 
N-body halo, we can study its reaction in the adiabatic 
limit, as well as for an abrupt change of the potential. In 
§2 we describe the model runs for different disk parame- 
ters and contraction times. In §3 we discuss the results 
and in §4 we draw our conclusions and present suggestions 
for future work. 

2. THE SIMULATIONS 

The N-body halo is set up according to the distribution 
function devised by Hernquist (1990, 1993). Its density 
distribution is p{r) = M;ir;,/[27rr(r -f VhY], where Mh is 
the total mass and is the scale length. The Hernquist 
halo has the same p on dependence in the center as the 
universal dark matter profile found by Navarro, Frenk, & 
White (1997), however with a finite total mass as the den- 



sity in the outermost regions decreases as p oc 



The 



baryonic component is represented by analytic external 
potentials of a disk and a central point mass. We used 
the potential of an exponentially thin disk according to 
Dehnen & Binney (1998). Its strength depends on the 
ratio of the disk mass Md to the disk scale length r^. Fol- 
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lowing Navarro, Eke, & Frcnk (1996), we started with a 
very large disk scale length and kept the disk mass con- 
stant throughout the simulation. The disk scale length is 
contracted linearly with time, where the contraction rate is 
a free parameter. Large contraction time scales compared 
to the dynamical time scales correspond to the adiabatic 
limit. In the limit of zero contraction time the halo will 
go through a phase of violent relaxation. The halo-disk 
system was allowed to relax after the contraction phase 
for several dynamical time scales. We applied a massive 
and a low-mass disk model with 20 % (MD) and 5 % (LD) 
of the total dark halo mass, and with a typical final scale 
length of 0.14 rh- 

In order to test the dependence on the concentration 
of the potential we used a softened point mass potential 
following 

^ ^ GMp arctan(r/rp) ^ 

Similar to the disk case we kept its mass Mp constant and 
shrank the smoothing length to a final value of rp — 0.03 
(test case PC). We also tested the case where Vp was kept 
constant and Mp was allowed to grow (PG). For all cal- 
culations we chose G = Mh = rh = 1, where G is the 
gravitational constant. Simulations were performed with 
halos represented by N = 1 x 10'', 8 x 10'', 2 x 10^ and 
1 X 10^ particles in order to test the dependence of the 
results on the numerical resolution. All model parameters 
are listed in Table 1. The gravitational softening length e 
was chosen according to the criterium of Mcrritt (1996). 
We used a time step at least a hundred times smaller than 
the dynamical time scale we would expect at the half mass 
radius for each model. The timestep, adequate for the runs 
with the highest number of particles, was not enlarged for 
lower particle number runs. We used a newly developed 
tree code WINE (Wetzstein et al., in preparation) in com- 
bination with special purpose hardware GRAPE-5 (Kawai 
et al. 2000) at the MPIA, Heidelberg. The refined force 
accuracy criterium of Salmon & Warren (1994) guarantees 
that the absolute force error stays below the precision of 
the GRAPE hardware. 

3. RESULTS 

The four top panels in Fig. 1 show the final density pro- 
files of the contracted halos for two resolutions for the mas- 
sive and the light disk, respectively (models LDl, LD2 and 
MDl, MD2). The error bars show the Poissonian error. 
The softening length is indicated with an arrow on each 
plot. In the innermost parts of the halo, i.e. r < 2e, the 
density is infiuenced by at least two effects: The softening 
length and fluctuations in the density due to small particle 
numbers. To show the effect of the fluctuations we plot the 
density profile of four consecutive dumps taken shortly af- 
ter the contraction has been completed. This is important 



as the fluctuations sometimes exceed the 1 a Poissonian 
error. At radii larger than two softening lengths the ana- 
lytical approximation of BW84 gives a very good account 
of the matter distribution of the halo for every disk model. 
The two panels at the bottom of Fig. 1 show the point 
mass case which is also in very good agreement with the 
theoretical profile, independent of the way we grow the ex- 
ternal potential. However, the density is somewhat lower 
than theoretically predicted inside the scale radius Tp. 

Fig. 2 illustrates for the LD case that the agreement 
between numerical simulation and theoretical profile im- 
proves with increasing particle number, as we are able to 
probe deeper into the center of the dark matter halo. 

In Fig. 3 we examine the response of the halo to an 
abrupt addition of the external potentials (models MDV, 
LDV and PV). The profile with the light disk (top left) 
still matches the theoretical curve. Even in the case of 
the massive disk (bottom left), the predicted density dis- 
tribution agrees well with the numerical model although 
we find a somewhat lower density inside a radius of one 
disk scale length. Adding the point mass potential instan- 
taneously (bottom right) leads to a density distribution 
that deviates strongly from the adiabatic prediction. The 
top right panel in Fig. 3 shows the time evolution of the 
1% mass shell radius in each case. For the LD models 
the mass shells stabilize at more or less the same radius, 
though the detailed path of contraction is very different 
for LD2 and LDV. For the more centrally concentrated 
potential (PV) this is not the case. It stabilizes at radii 
further outside than in the adiabatic case. Consequently 
these density profiles are less concentrated in the center 
than in the adiabatic case. 



4. CONCLUSIONS 

We find that the adiabatic approximation is very robust 
in the parameter space that is occupied by normal spiral 
or dwarf galaxies. This remains valid even in cases where 
the contraction timescales are short compared with the 
dynamical time scales. Only for a point mass potential, 
a significant difference occurs. We conclude that the va- 
lidity of the adiabatic approximation depends mainly on 
the concentration of the added potential. A better under- 
standing of the underlying physics will require a detailed 
analysis of the orbital distribution of the particles prior 
and after the contraction phase. We will investigate this 
question in a subsequent paper. 

A program to calculate disk potentials was kindly pro- 
vided by Walter Dehnen. We want to thank Hans- Walter 
Rix, Matthias Steinmetz and Chris Gottbrath for helpful 
comments and discussions. 
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Table 1 

Model nrxs 



Model 


Mass/scale length 


N 


e 


dt [10-3 


] type 


LDl 


0.05/0.14 


1 X 10^ 


0.03 


6.25 


adiabatic 


LD2 


0.05/0.14 


2 X 105 


0.007 


6.25 


adiabatic 


LD3 


0.05/0.14 


1 X 10^ 


0.0035 


6.25 


adiabatic 


LDV 


0.05/0.14 


2 X 10^ 


0.007 


6.25 


violent 


MDl 


0.2/0.14 


1 X 10" 


0.03 


1.57 


adiabatic 


MD2 


0.2/0.14 


2 X 10^ 


0.007 


1.57 


adiabatic 


MDV 


0.2/0.14 


2 X 10-5 


0.007 


1.57 


violent 


PG 


0.2/0.03 


8 X 10-* 


0.01 


0.78 


point grow 


PC 


0.2/0.03 


8 X 10^ 


0.01 


0.78 


point contract 


PV 


0.2/0.03 


8 X 10^ 


0.01 


0.78 


point violent 


Note. 


— Abbreviations for the model 


runs are LD = 


Light Disk, MD = 



Massive Disk, P = Point Mass, V = Violent, G = Mass Growth, C = Softening 
Length Contraction. Model parameters are N = number of particles and dt = 
integration time step. Second column shows mass to scale length ratios for point 
mass and disk models. 




Fig. 1. — Final density profiles of contracted dark halos. The upper two panels show the results for the models LDl and LD2, the following 
two for the models MDl and MD2. The last two panels show the results for the point mass. In the left panel the contracting scale length 
case and in the right panel the growing mass case is shown. Arrows indicate the used softening length. 
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Fig. 2. — Final density profiles of LD1,LD2 and LD3 are shown to examine the effects of numerical resolution. 
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Fig. 3. — Final density profiles after abrupt addition of external potential. LDV top left, MDV bottom left and PV bottom right. In the 
top right panel a comparison of the evolution of the 1% mass shell radii between the adiabatic and the violent case of the models LD and P 
is shown. 



